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Abstract 

The class of stochastic Runge-Kutta methods for stochastic differential equations 
due to RoBler is considered. Coefficient families of diagonally drift-implicit stochas- 
tic Runge-Kutta (DDISRK) methods of weak order one and two are calculated. 
Their asymptotic stability as well as mean-square stability (MS-stability) proper- 
ties are studied for a linear stochastic test equation with multiplicative noise. The 
stability functions for the DDISRK methods are determined and their domains of 
stability are compared to the corresponding domain of stability of the considered 
test equation. Stability regions are presented for various coefficients of the fami- 
lies of DDISRK methods in order to determine step size restrictions such that the 
numerical approximation reproduces the characteristics of the solution process. 
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1 Introduction 



Numerical methods are an important tool for the calculation of approximate 
solutions of stochastic differential equations (SDEs) which possess no analyti- 
cal solution formula. Therefore, many approximation schemes have been devel- 
oped in recent years and much research has been carried out to develop deriva- 
tive free stochastic Runge-Kutta (SRK) type methods [Z|17jl8|19|Z7] . Similar 
to the well understood deterministic setting of ordinary differential equations 
(ODEs), one has to pay much attention to the stability properties of the solu- 
tion as well as of the numerical approximations. Therefore, implicit methods 
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have been proposed for the strong pathwise approximation of solutions of SDEs 
and their stability has been analyzed [2.9, l6|25] . However, for the approxima- 
tion of moments of the solution process special numerical methods converging 
in the weak sense have to be applied (see, e.g., [TT(12|ll3|ll41ll7|ll8|ll9f2Uf21f27] ) 
and a stability analysis has to be carried out similar to that for strong ap- 
proximations [9|13|26j . In the present paper, we present families of first and 
second order diagonally drift-implicit SRK (DDISRK) methods for the weak 
approximation of SDEs contained in the class of SRK methods proposed by 
Rdfiler [20] ■ Further, we analyze their asymptotic stability and mean-square 
stability for linear test equations with multiplicative noise. Finally, the regions 
of stability of the DDISRK methods are compared to the regions of stability 
of the linear test equation. Thus, in Section [2] we consider the class of SRK 
methods and coefficients families for weak order one and order two DDISRK 
methods are presented. Then, we discuss the concepts of stability for solutions 
of SDEs and for numerical approximations in Section [3] and Section HJ respec- 
tively. In Section |5l some numerical experiments are carried out in order to 
justify our theoretical results. 

Let (Q, J 7 , P) be a probability space with a filtration (J"i) t > which fulfills 
the usual conditions and X = [to,T] for some < t < T < oo. Then, let 
X = (X t ) te x denote the solution process of an Ito SDE 

dX t = a{t, X t ) dt + b{t, X t ) dW t , X t0 = x , (1) 

for t G X where a : X x R d -»■ R d is the drift and b : X x R d -»■ R dxm is the 
diffusion, (W t )t>o is a m-dimensional Wiener process and a Xi -measurable 
initial condition x independent of W t —W to for t > t such that E(||xo|| 2r ) < oo 
for some r G N. The jth column of the d x m-diffusion matrix b = (b 1 ^) will 
be denoted by V in the following. Further, we suppose that the conditions of 
the existence and uniqueness theorem [TT] are fulfilled for SDE (pQ). 

In the following, we consider time discrete approximations Y h = (Y t ) t£ x h w.r.t. 
a constant step size h = for some N E N and Ih = {to, h, . . . , t^} 

where t n = t + n h for < n < N . As usual, we also write Y n = Y tn for 
< n < N. Further, let C^(M d ,M) denote the space of all g G C"(M d ,M) 
fulfilling a polynomial growth condition [11] . 

Definition 1.1 A time discrete approximation Y h converges weakly with or- 
der p > to X as h -»■ at time t G X h if for each f G Cp (p+1) (R d , R) exists 
a constant Cf, which does not depend on h, and a finite 5q > such that for 
each h G ]0, 5q[ 

\E{f{X t ))-E(f{Y t ))\<C f hP. (2) 
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2 Diagonally Drift-Implicit Stochastic Runge— Kutta Methods 



For the weak approximation of the solution of the Ito SDE (pQ), we consider 
the class of SRK methods introduced by Rofiler [20] • Then, the d-dimensional 
approximation process Y h with Y n = Y tn for t n G is given by the following 
SRK method of s- stages with Y = x and 

s 

Y n+1 = Y n + ^2 a i a (tn + c[ 0) h n , H> 0) ) h n 

i=l 

s m f 

+ E E (A W W + A (2) -W) b k (t n + c?h n , /ff ) (3) 



+ e e (a (3) h),n + vM b\ tn + c^h n , m 

i=l k=l 

for n = 0, 1, . . . , iV — 1 with stage values 

HP>= Y n + £ 4 0) a(t n + cf Hf) /i n 



+ EE s S ,)6i (*«+4 1)/i "^f)W 
i=i ;=i 

H\ k) = Y n + j2 4j <tn + cfh n , iff) h n 

3=1 



3=1 

H\ k) = Y n + J2 4f a(t n + cf>hn, iff) h n 

3=1 

s m f 

+ T,T,Bifb\t n + cfh ni H^ 1 ^ 



7=1 1=1 V IL n 

l^k 

for % = 1, . . . , s and k = 1, . . . , m. The random variables of the method are 
defined by 

{\{I(k),nI{f),n ~ Vhnl(k),n) H k < I 

l(I(k),rJ(l),n + VK,I(l),n) H I < k (4) 

Whn-K) if fc = Z 

for 1 < k, I < m with independent random variables I(k),n for 1 < k < m and 
I{k),n for 1 < k < m—1 and < n < N. Thus, only 2m— 1 independent random 
variables have to be simulated for each step. In the following, we choose I<k), n 
as a three point distributed random variable with P(I^),n = ±a/3 h n ) = | 
and P(J(fe) jn = 0) = |. The random variables I(k),n are defined by a two point 
distribution with P(I( k ^ n = ±y/h) = |. 
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The main advantage of this class of SRK methods is the significant reduction 
of complexity compared to present SRK methods in recent literature, because 
the number of stages does not depend on the dimension m of the driving 
Wiener process [20J. We denote by a = (a^ ) and = ($ k) ) for 1 < k < 4 
the corresponding vectors of weights and by = {A^) and B^ = (B^) 
for k — 0,1,2 the corresponding coefficients matrices. Then, the coefficients 
of the SRK method ([3]) can be represented by an extended Butcher array: 



c (0) 


A <P) 


B(°) 








B(i) 






A m 


B (2) 






a T 















Weak order one and two conditions for the SRK method ([3]) have been cal- 
culated in [2U] by applying the colored rooted tree theory due to RoBler in- 
troduced in [117]. Now, let ps = P denote the order of convergence of the SRK 
method ([3]) if it is applied to an SDE and let pr> with pr> > ps denote the 
order of convergence if it is applied to a deterministic ODE, i.e., SDE ([Q) with 
6 = and we also write (pd,Ps) [HUB] - Since we are interested in SRK meth- 
ods which inherit good stability properties, we consider families of DDISRK 
methods which are diagonally implicit in the deterministic part of the scheme. 



2.1 Weak Order One DDISRK Methods 



Firstly, we consider weak order one DDISRK methods ([3]) with s = 1 stage |20j . 
However, in order to cover the stochastic theta method PTTT] , we also consider 
the case that the stage number is s = 2 for the drift function only, whereas it 
is still one for the diffusion function. Then, from the order conditions [20] it 
follows that the family of weak order one DDISRK methods is characterized 
by the Butcher table ([5]) with some coefficients c±, c 2 , c 3 , c 4 , c 5 G E. As an 
example, in the case of s — 1 stage we obtain for c\ = c 2 = c 3 = c 4 = c 5 = 
the explicit Euler-Maruyama scheme of order (1, 1) [TT]. For c\ = | and c 2 = 
C3 = C4 = C5 = we obtain the SRK scheme DDIRDI1 with s = 1 stage of 
order (2, 1), which reduces to the midpoint rule if it is applied to an ODE [5]. 
If we consider the case of s = 2 stages, then we get for c\ = 0, c 2 = 1 — 0, 
c 3 = e, c 4 = 1 and c 5 = 9 for some 9 E [0,1] the SRK scheme DDIRDI2 
of order (2, 1) which coincides with the stochastic theta method Plllf23f24j . 
Further, for s = 2 stages with c\ = c 3 = 3+ 6 , c 2 = — c 5 = | and c 4 G K 
we get DDISRK schemes of order (3, 1) which are A-stable in case of ODEs 
[5]. Especially, in the case of c 4 = | we denote the scheme as DDIRDI3 in the 
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following. Note that the schemes DDIRDI1 and DDIRDI2 for 6 — | are also 
A-stable if they are applied to ODEs [5]. 

2.2 Weak Order Two DDISRK Methods 

Next, we consider weak second order DDISRK methods ([3]) with s = 3 stages. 
Here, we claim that a% = and — for 1 < i,j < 3 in order to reduce 
the computational effort. Since in this case the third stage does not 
matter anymore, we let A$ = B$ = for 1 < j < 3. Then, we can obtain 
oil — OL2 = | from some order conditions of weak order two [20] and from 
the classification given in [1] in the case of an explicit SRK scheme. On the 
other hand, since we assume = 0, all conditions for A^ to satisfy are 
a 7 A^e = | only. Therefore, we can consider arbitrary coefficients A[j as 
long as a T A^e = | is fulfilled for o>i = «2 — \- As a result of this, the weak 
order two DDISRK schemes ([3]) are given by the infinite coefficients family 
(EJ) with ci,C2 G 1 and C3, C4 G K \ {0}. Clearly, one has to solve 2 (in general 
nonlinear) systems of equations from the stage values, each of dimension d, for 
the DDISRK method if c\ 7^ and C2 7^ 0. Therefore, as in the deterministic 
setting, some simplified Newton iterations have to be performed in each step 
in order to solve the nonlinear system of equations [3]l5] . As an example, for 
Ci = c 2 = 3+ /^ and for all 03,04 6 R \ {0} we obtain an DDISRK scheme of 
order (3,2) which is A-stable if it is applied to a deterministic ODE [5]. 



3 Stability analysis for SDEs 

For SDEs several stochastic stability concepts have been proposed in literature, 
see e.g., [9|10fllfl3|rT5|2 2 23 24 26j and the literature therein. In the following, 
we consider SDE (JTJ) with a steady solution X t = such that a(t, 0) = b(t, 0) = 
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holds, which is also called an equilibrium position. Suppose that there exists 
a unique solution X t = X(t;t ,x Q ) for all t > t and for each nonrandom 
initial value Xq under consideration. Then, stochastic stability can be defined 
as the stochastic counterparts of stability, asymptotic stability and asymptotic 
stability in the large for ODEs [T|5fTTj . 

Definition 3.1 Let X be the solution of the scalar ltd SDE (T7J). Then, the 
equilibrium position of the SDE is said to be 

(i) stochastically stable if for all e > and t > holds 



lim P sup \X(t;t ,x )\ > e =0, 
*o->0 \t>to J 

(ii) stochastically asymptotically stable if (TJj holds and if 

lim p( lim \X(t;t ,x )\=0) = 1, 

(Hi) or stochastically asymptotically stable in the large if (TJj holds and if for 
all x G M holds 



P[\im\X(t;t ,Xo)\ = 0) =1. 



Further, stability analysis involving the pth moments of the solution process 
is also widely considered, see e.g. |T|[9]llO|llHfT3] . 

Definition 3.2 Let X be the solution of the scalar ltd SDE (T7J). Then, the 
equilibrium position of the SDE is said to be 
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(i) stable in the pth-mean if for every e > and to > there exists a 
5 = 5(to, e) > such that for all t > t and \xq\ < 5 

E(\X(P,to,x )\*)<e, 

(ii) asymptotically stable in the pth-mean if ||) holds and if there exists a 
= ^o(^o) > such that for all \xq\ < So 

\imE(\X(t;to,x )\ p ) = 0- 

t— >oo 

The most frequently used cases in applications are p = 1 and p — 2, i.e., 
stability in mean (M-stability) and mean-square stability (MS-stability). In 
the present paper, we will focus on asymptotic stability in the large and MS- 
stability for a linear test equation with multiplicative noise [8f9~|10|l23] 

dX t = XX t dt + pX t dW t (7) 

for t > t and with some constants A, \x G C and with a nonrandom initial 
condition X to = xq G M\{0}, which reproduces the dynamics of more complex 
SDEs better than in the case of additive noise [71TTT] . The exact solution of ([7j) 



can be calculated as X t = xq exp((A — ^p )(t — t ) + p(W t — W to )) which is 



stochastically asymptotically stable in the large [23] if 

lim \X t \ =0 with probability 1 Sft(A - \p 2 ) < . 



t— >oo 



We calculate that \X t \ p = \x \ p exp(p3^(A - \p 2 ){t - t ) + p$t(p)(W t - W to )) 
which yields E(|X 4 | P ) = |x | p exp(pK(A - \p 2 ){t - to) + \p 2 {^{p)) 2 {t - to)). 
Then the pth-mean stability domain where SDE (J7j) possesses an equilibrium 
position can be determined as follows: 

limE(|X t n = <=> 2 5R(A)-% 2 )+p(%)) 2 <0. (9) 

Thus, the equilibrium position of SDE (J7J) is asymptotically MS-stable if 

lim E(\X t \ 2 ) = ^ 2$l(X) + \p\ 2 <0 (10) 

t— >OQ 

for the coefficients X, p G C (see, e.g., [9|23f26] ). We remark that due to 
9ft(2A— p 2 ) < 2 K(A) + |/i| 2 MS-stability always induces asymptotically stability 
in the large. Further, for p = the stability condition (Q reduces to the well 
known deterministic stability condition 3?(A) < 0. 



4 Numerical stability of SRK methods 



We are now looking for conditions such that a numerical method applied to 
SDE (J7J) yields numerically stable solutions. Therefore, we say that the method 
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is numerically asymptotically stable or MS-stable if the numerical solutions 
Y n satisfy lim^oo \Y n \ = with probability one or lim^oo E (|F„| 2 ) = 0, 
respectively. If we apply the numerical method to the linear test equation 
flTJ, then we obtain with the parametrization h = Xh and k = /jy/h PIT3"] a 
one-step difference equation of the form 

n 

Y n+1 = R n (h, k) Y n = J] lUih, k) Y . (11) 

i=0 

with a stability function R n (h,k). The domain of asymptotic stability of a 
numerical method can be determined by the following lemma [9]: 

Lemma 4.1 Given a sequence of real-valued, non-negative, independent and 
identically distributed random variables (\R n (h, fc)|) n6 N ? consider the sequence 
of random variables (\Y n \) ne ^ defined by (Tl\) where \Y \ ^ with probability 
1. Suppose that the random variables log(R n (h, k)) are square-integrable. Then 

lim \Y n \ = , with probability 1 <^> E(\og(R n (h, k))) < . (12) 



We call the set H AS = {(h,k) G C 2 : E(\og(R n (h,k))) < 0} C C 2 the do- 
main of asymptotical stability of the method. Note that one can also find 

2 

some alternative parameterizations like k = — y- in the literature [T|23f26] . 
Analogously, if we calculate the mean-square norm z n = E(|l^| 2 ) then we 
obtain a one-step difference equation of the form z n+1 = R(h, k) z n where 
R(h,k) = E(\R n (h, k)\ 2 ) is called the MS-stability function of the numerical 
method. Thus, we obviously yield MS-stability, i.e. z n — > as n — > oo, if 
R(h,k) < 1. The set Kms = {{h,k) E C 2 : R(h, k) < 1} C C 2 is called the 
domain of MS-stability of the method. 

Especially, the domain is called region of stability in the case of (h, k) e M. 2 
|23j . The numerical method is said to be ^-stable if the domain of stability 
of the test equation ([7]) is a subset of the domain of numerical stability. Since 
the domain of stability for A, /i e C is not easy to visualize, we restrict our 
attention to figures presenting the region of stability with A, fi G R in the h-k 2 
plane. Then, for fixed values of A and /i, the set {(A h, // 2 h) C M 2 : h > 0} is a 
straight ray starting at the origin and going through the point (A, /j, 2 ). Clearly, 
varying the step size h corresponds to moving along this ray. For A, [i G M, 
the region of asymptotical stability for SDE ([7]) reduces to the area of the 
h-k 2 plane with the /i-axis as the lower bound if h < and with k 2 > 2h 
as the lower bound if h > whereas the region of MS-stability for SDE (J7J) 
reduces to the area of the h-k 2 plane with the /i-axis as the lower bound and 
k 2 < —2h as the upper bound for h < 0. 

In the following, all figures presenting regions of stability for some numerical 
method under consideration are plotted by the software Mathematica. The 
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regions of numerical asymptotically stability and MS-stability are indicated 
by two dark-grey tones whereas the regions of MS-stability are more dark 
than the regions of asymptotical stability Further, the corresponding regions 
of stability for the test equation ([7]) are filled by two light-grey tones whereas 
again the regions of MS-stability are more dark than the regions of asymptot- 
ical stability. In all presented figures, the regions of MS-stability are a subset 
of the regions of asymptotical stability 

4-1 Stability of Order One DDISRK Schemes 

We consider the family of order one DDISRK schemes §5§ with coefficients 
ci, . . . , C5 GR. If we apply these schemes to the linear test equation ((Tj) then 
we obtain the difference equation 

Y n+l = Y n + (l-c 5 )\h H{ 0] + c 5 Xh Hi 0) + fi I (1) , n H{ 1} (13) 

with the stage values 

H[ 0) =Y n + Cl \hH[ 0) 

Hf ] =Y n + c 2 Xh h[ 0) + c 3 Xh #< 0) + c 4 n I (1)) „ (14) 
H? = Y n 

where the implicit equations for h[ ^ and H% can be solved in the case of 
1 — Ci A h 7^ and 1 — C3 A h 7^ 0, which is fulfilled for step sizes h 7^ ^ if 

ci 7^ and h 7^ ^ if c 3 7^ 0. With h = A h and k = fj,y/h let 

r _ 1 | h - c 3 h 2 + c 5 (C2 + c 3 - ci) fr 2 S - ° 4 ° 5 ^ fc I 

(1 — Ci /l)(l — C3 A) 1 — C 3 /l 

Then, we can write ffl3|) by the recursion formula Y n+ i = R n (h, k) Y n with 
the stability function R n (h\ k) = V + hr 1 ! 2 E /(i), n for n = 0, . . . , JV — 1. Since 
the SRK schemes (jHJ) are of weak order one, we can substitute the tree point 
distributed random variables I(j), n by two point distributed random variables 
i(j), n for 1 < j < m in (EJ) and consider R n (h, k) = T + /i -1 / 2 E instead. 
Now, we analyse the asymptotic stability of the SRK schemes dSJ) by apply- 
ing Lemma 14.11 Further, in order to analyse the MS-stability, we calculate 
the mean-square norm z n = E(|Y„| 2 ). Then, we obtain the recursion formula 
z n+ i = R(h, k) z n with the MS-stability function R(h, k) = \T\ 2 + |E| 2 . 

Proposition 4.2 For SDE |?|] with A, /x G C, the SRK schemes (EJj are 

(i) numerically asymptotical stable if \T 2 — 3 E 2 | |T| 4 < 1 and in the case 
that the random variables n are replaced by Iu\ n for 1 < j < m, if 

|r 2 -s 2 |<i, 
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Fig. 1. Asymptotical and mean-square stability region for DDIRDI1 with two point 
and three point distributed random variables in the left and right figure, respectively. 



k"2 k A 2 




Fig. 2. Asymptotical and mean-square stability region for DDIRDI2 with two point 
and three point distributed random variables in the left and right figure, respectively. 

(ii) numerically MS-stable i/|r| 2 + |£| 2 <l. 

Here, we have to point out that the distribution of the random variables 
used for the numerical method has significant influence on the domain of 

asymptotical stability. As an example, for DDIRDI1 we have F = | + i^ , £ = k 

and we calculate K AS = {(h,k) G C 2 : |(1 + \h) 2 - (1 - \h) 2 k 2 \ < \l-\h\ 2 } if 
the random variables I(j), n are used and IZas = {(h, k) G C 2 : |(1 — \h) 2 (\ + 
\h) 2 (l — 3A; 2 )||1 + |/i| 4 < |1 — ^h\ 4 } if I(j),n are used. In both cases, we get 
n MS = {(h,k) G C 2 : 2^(h) + ((1 - \%l{h)) 2 + \{%h)) 2 ) \k\ 2 < 0}. The 
corresponding regions of stability are presented for both cases in Figure [TJ For 
DDIRDI2 we get Y = 1+< ; l ~f l and E = -^k + k, 9 G [0, 1]. Then, for 9 = ± 

1— oh L—Bh ^ 

follows Has = {(h,k) G C 2 : |(1 + \h) 2 - k 2 \ < |1 - \h\ 2 } if %„ are used, 
^5 = {(h,k) G C 2 : |(l-|/i) 2 ((l + i/ i ) 2 -3A: 2 )||l + i/i| 4 < |l-|/i| 4 } if / ov are 
used and K M s = {{h, k) G C 2 : 291(h) + \k\ 2 < 0}. Thus, the scheme DDIRDI2 
with 9 — | is A-stable w.r.t. MS-stability and the corresponding regions 
are presented in Figure [2] (see also [91123] ). Analogously, we can calculate the 
domains of stability for DDIRDI3 which are presented in Figure |3j For all 
considered schemes, we can see the influence of the random variables used by 
the scheme to the domain of asymptotical stability. 
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Fig. 3. Asymptotical and mean-square stability region for DDIRDI3 with two point 
and three point distributed random variables in the left and right figure, respectively. 

4-2 Stability of Order Two DDISRK Schemes 



Next, we apply the DDISRK method ([3]) with the coefficients (EJ) to the linear 
test equation (jTJ). Then we obtain the difference equation 

Y n+1 = Y n + ±\hH[ 0) + ±\hHi 0) 

+ ( 1 - ^) " f « H ® + h * ha H P + hi) (is) 



2c 3 P Vh 2 2c 3 P Vh 3 



with stage values 

H^ = Y n + Cl XhH^ 

= Y n + {l-c l -c 2 )Xh H[ 0) + c 2 Xh if 2 (0) + n J (1) 

= Y n (16) 

= Y n + c\ A /i H[ 0) + c 3 fiVh H[ 1] 
HP = Y n + c\Xh H[ 0) -c 3f iVh H[ 1] 

where the values do not appear due to m — 1 and = 0. Suppose that 
1 — c\ A h 7^ and that 1 — c 2 A /i ^ which can always be fulfilled for step 
sizes h with h ^ ^ and ft, 7^ Then the implicit equations for and 

can always be solved. Let with h = Xh and k = \x\fh 

r = i + U( ^: Cl - C2)A2 , e = A = V 

(l- Cl /i)(l-c 2 /i) (l- Cl /i)(l-c 2 /i) 2 

Then, we yield for (JT5|) the recursion formula Y^ +1 = R n (h,k)Y n with the 
stability function R n (h, k) = T — A + /zr 1 / 2 £ i(i), n + ^ 1 A n . In order to 
analyse the asymptotic stability of the weak order two DDISRK schemes ([6]) we 
apply again Lemma 14.11 For the determination of the domain of MS-stability, 
we calculate the mean-square norm of (TToT) . Then, we obatin the recursion 
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Fig. 4. Asymptotical and mean-square stability regions for RI6 with c\ = C2 = on 
the left and for DDIRDI4 with c\ = and C2 = \ on the right. 

z n+ i = R(h,k)z n with MS-stability function R(h, k) = \T\ 2 + |£| 2 + 2|A| 2 . 
Both stability functions i? n (/i,fc) and R(h,k) depend only on the coefficients 
Ci and C2 of the scheme, i.e. the coefficients C3 and C4 are not relevant for the 
stability in the case of the scalar linear test equation (J7J). 

Proposition 4.3 For SDE |7|) with X, fx G C, the SRK schemes (T2|) are 

(i) numerically asymptotical stable if \{T + 2 A) 2 — 3S 2 | |T — A| 4 < 1, 

(ii) numerically MS-stable if \T\ 2 + |S| 2 + 2|A| 2 < I. 

If we choose c\ = C2 = and C3 = C4 = 1 in (jSJ), then we yield the explicit SRK 
scheme RI6 calculated in [20] which coincides for m = 1 with the SRK scheme 
due to Platen [ll.18.26j. If we apply Proposition 14.31 for RI6, then we obtain 
n AS = {(h,k) G C 2 : \(l + h+fh 2 + k 2 ) 2 -3(l + h) 2 k 2 \\l + h+fh 2 -lk 2 \ 4 < 1} 

and 11ms = {(h, k) G C 2 : |1 + h + \h 2 \ 2 + |1 + h\ 2 \k\ 2 + ||A;| 4 < 1}. The 
corresponding regions of stability are given in Figure HI Further, we can choose 
c\ = 0, C2 = \ and e.g. C3 = C4 = 1 which defines the scheme DDIRDI4. Then, 
the DDIRDI4 scheme ([3]) is an advancement of the stochastic theta method 
DDIRDI2 with 9 = \. However, for DDIRDI4 we get 1Z AS = {(h,k) G C 2 : 

K 1 + lA I + fc2 ) 2 -3(l + ^) 2 fc 2 ||l + I A I -|A; 2 | 4 < 1} andTZ MS = {(h,k) G 
^ 2 : |l + 7zbrl 2 + |l + T^E I 2 |^| 2 + |l^l 4 < 1}- The regions of stability are given 

1 2 h 1 2 

in Figure HI Here, we can see that the good stability properties of the order one 
scheme DDIRDI2 are not carried over to the second order scheme DDIRDI4. 
Therefore, we are looking for further second order DDISRK methods with 
some better stability qualities. 

It is usual to consider singly diagonally implicit Runge-Kutta methods for 
ODEs where all coefficients A$ are equal. Therefore, we assume that c\ = C2 
for the schemes ([6]) in the following. Then, the domains of stability are TZas — 

{(*,*) 60:|(l + ^^+t>)'-3(l + ^)»f||l + 4f^-l*>|«< 

1} and Tl MS = {(h,k) € C 2 : |1 + I 2 + |1 + |^I 2 I*I 2 + ^l*l 4 < 

1} which depend on the coefficient c\ of the scheme. In the following, we 
consider various values for the parameter c\ of the DDISRK scheme (jSJ) and 
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Fig. 5. Stability regions for DDIRDI5 with c\ = \ on the left and c\ = \ on the 
right. 
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Fig. 6. Stability regions for DDIRDI5 with c\ = ^ + on the left and c\ = 1 on 
the right. 
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Fig. 7. Stability regions for DDIRDI5 with c\ = § on the left and c\ = 3 on the 
right. 

we analyze the stability domain for A, /i G R. Therefore, we choose c\ G 
{|, |, | + tjV^, 1, §, 3}. Especially, we consider the case of C\ = c 2 = \ + |\/3 
and C3 = C4 = 1 which we denote as the scheme DDIRDI5. Then, and a 
coincide with the coefficients of the well known deterministic SDIRK scheme 
which is A-stable and attains order po = 3 for deterministic ODEs [5]. The 
corresponding stability regions are presented in Figures [5H3 



5 Numerical Experiments 



We compare the efficiency of the proposed second order DDISRK schemes 
DDIRDI4 and DDIRDI5 with the second order drift-implicit SRK scheme 
DIPL1WM due to Platen ([H|, p. 501). Therefore, we take the number of 
evaluations of the drift function a, of the diffusion functions fr 7 ', j — 1, . . . , m, 
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Fig. 8. Error vs computational effort with double logarithmic scales for SDE (|17|) 
on the left and SDE (fTHj) on the right hand side. 



and the number of random numbers needed each step as a measure of the 
computational complexity for each considered scheme. As a first example, we 
consider for d — m — 1 the Ito SDE 



dX t = [\X t + sjxj+l) dt + yJx? + ldWt , X 



0, 



(17) 



on the time interval / = [0, 2] with solution X t = smh(t + W t ). Here, we choose 
f(x) = p(arsinh(x)), where p(z) = z 3 — 6z 2 + 8z is a polynomial. Then we 
calculate that E(f(X t )) = t 3 — 3t 2 + It which is approximated at time t = 2 
with step sizes 2 _1 ,...,2~ 4 and 10 8 simulated trajectories. The results are 
presented on the left hand side of Fig. [81 As a second example, we consider a 
nonlinear SDE with a 10-dimensional driving Wiener process 

dX t = Xt dt + ^x t + l -dwl + ^X t + \dW? + l^ + Idw; 



10 



15 , 

— \ X t + — dWf + —\ X t + - dWf + —\ X t + - dW' t 
40 V 20 * 25 V 2 * 20 V 4 * 



+ — « Ax* + - dWf + — Jx t + — dW. 9 + — Jx t + — diy, 10 , X = 1. 
15 V 5 * 20 V 10 4 25 V 20 * 



4 
1 

25 
1 

25 



20 

T 

2 



5 
1 

20 



(18) 



with non-commutative noise. Here, we approximate the second moment of 



68013 



+ 



68013 



+ 1) exp 



731453 



the Solution E(X t ) — 1462 9060 1 V14629060 1 "^ wv ^^360000' 

simulated trajectories with step sizes 2°, . . . , 2~ 3 . The results are presented 



t) at time t — 1 by 



10 b 



on the right hand side of Fig. [HI Here, the schemes DDIDRI4 and DDIDRI5 
perform impressively better than the drift-implicit scheme DIPL1WM [TTj . 
This is a result of the reduced complexity for the new class of efficient SRK 
schemes due to Rofiler [20] which becomes significant especially for SDEs with 
high-dimensional driving Wiener processes. 



Next, we verify the theoretical results for the domains of stability of the pro- 
posed SRK methods by numerical experiments. Therefore, we consider the test 
equation (JTj) with parameters A = —200, fi = Vo and with initial value X = 1 
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Fig. 9. Asymptotical and MS-stability analysis for RI6, DDIRDI1 and DDIRDI5. 

on the time interval [0, 10]. We apply the second order explicit SRK scheme 
RI6 [20J, the order one DDISRK scheme DDIRDI1 with two point as well 
as with three point distributed random variables and the order two DDISRK 
scheme DDIRDI5. We denote by DDIRDI1-2P the scheme DDIRDI1 if 2 point 
distributed random variables -%), n are used instead of im, n - 

In order to analyse the numerical asymptotically stability, a single approxima- 
tion trajectory is simulated with each scheme under consideration for the step 
sizes h = 0.01, h = 0.1, h = 0.2 and h = 0.5. Then, we obtain the following 
theoretical results due to Proposition 14.21 and Proposition 14.31 the scheme RI6 
is asymptotical stable for the step size h = 0.01 and it is unstable for h = 0.1, 
h = 0.2 and h = 0.5. DDIRDI1 and DDIRDI1-2P are stable for h = 0.01 and 
h = 0.1. In the case of h = 0.2 only DDIRDI1-2P is stable while DDIRDI1 
is unstable. Further, DDIRDI1 and DDIRDI1-2P are unstable for h = 0.5. 
Finally, the scheme DDIRDI5 is asymptotical stable for h = 0.01, h — 0.1 and 
even for h = 0.2, however it is unstable for h = 0.5. The numerical results 
for a single trajectory \Y n \ are plotted with logarithmic scale to the base 10 
versus the time on the left hand side of Fig. [91 We remark that the results for 
DDIRDI1 with step size h = 0.01 tend to zero after two steps and are thus 
not visible in Fig. [9j 

For the analysis of the numerical MS-stability, the value E(|X t | 2 ) is approx- 
imated by Monte Carlo simulation with 10 4 independent trajectories for the 
step sizes h = 0.005, h = 0.01, h = 0.1 and h = 0.5. Proposition 14.21 and 
Proposition 14.31 give the following results: RI6 is MS-stable for h = 0.005 and 
MS-unstable for all other considered step sizes. DDIRDI1 and DDIRDI1-2P 
are MS-stable for h = 0.005 and h = 0.01, however MS-unstable for h = 0.1 
and h = 0.5. Further, DDIRDI5 is MS-stable for step sizes h = 0.005, h = 0.01 
and even for h = 0.1 and MS-unstable for h = 0.5. The corresponding nu- 
merical results of E(|Y^| 2 ) are presented with logarithmic scale to the base 10 
versus the time on the right hand side of Fig. [9j Again, the numerical results 
exactly confirm our theoretical findings for the domains of stability. 
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